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Retinoic acid (RA) triggers physiological processes by activating heterodimeric transcription factors 
(TFs) comprising retinoic acid receptor (RARa, p, y) and retinoid X receptor (RXRa, p, y) . How a 
single signal induces highly complex temporally controlled networks that ultimately orchestrate 
physiological processes is unclear. Using an RA-inducible differentiation model, we defined the 
temporal changes in the genome- wide binding patterns of RARy and RXRa and correlated them with 
transcription regulation. Unexpectedly, both receptors displayed a highly dynamic binding, with 
different RXRa heterodimers targeting identical loci. Comparison of RARy and RXRa co-binding at 
RA-regulated genes identified putative RXRa-RARy target genes that were validated with subtype- 
selective agonists. Gene-regulatory decisions during differentiation were inferred from TF-target 
gene information and temporal gene expression. This analysis revealed six distinct co-expression 
paths of which RXRa-RARy is associated with transcription activation, while Sox2 and Egrl were 
predicted to regulate repression. Finally, RXRa-RARy regulatory networks were reconstructed 
through integration of functional co-citations. Our analysis provides a dynamic view of RA 
signalling during cell differentiation, reveals RAR heterodimer dynamics and promiscuity, and 
predicts decisions that diversify the RA signal into distinct gene-regulatory programs. 
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Introduction 

Retinoic acid receptors (RARs) and retinoid X receptors (RXRs) 
are members of the nuclear receptor (NR) gene family of 
ligand-regulated transcription factors (TFs). RARs and RXRs 
form heterodimers that act as master regulators for multiple 
physiological processes, including embryogenesis, organo- 
genesis, immune functions, reproduction, and organ homeo- 
stasis (Mark et al, 2006). Apart from their impact on 
physiology, RARs and RXRs have major promise for therapy 
and prevention of cancer and other diseases, and several 
therapeutic paradigms have been established (Altucci et al, 
2007; Liby et al, 2007; Shankaranarayanan et al, 2009; 
de The and Chen, 2010; Zhang et al, 2010). 

The biological importance of the retinoid signalling system 
and its cancer therapeutic potential has inspired intense 
research that provided detailed insight in the structural 
basis of, and molecular events at the early steps of retinoid 
action. Mechanistically, the binding of a ligand facilitates 
the exchange between corepressor (CoR) and co-activator 
(CoA) complexes by allosterically altering receptor surfaces 
involved in these interactions. The recruitment of such 
epigenetically active and/or chromatin modifying complexes 



leads to chromatin structure alterations and post-translational 
modifications that ultimately regulate cognate gene programs 
(Gronemeyer et al, 2004; Rosenfeld et al, 2006). 

The retinoid signalling system is highly complex, as it 
comprises three RXR (RARa, (3 and y) and three RAR (RARa, 
(3 and y) subtypes expressed from distinct genes as multiple 
isoforms which act as heterodimers; in addition, RXRs can 
form heterodimers with a plethora of other NRs (Laudet and 
Gronemeyer, 2002) . While insight into (some of) the physio- 
logical functions of the various RAR and RXR subtypes has 
been obtained by exploiting mouse genetics (Mark et al, 
2006) our understanding of the cell physiological functions 
of these various subtypes is rather limited. The generation 
of subtype-selective ligands has provided important tools 
(de Lera et al, 2007), while the study of RAR subtype-deficient 
F9 embryonal carcinoma (EC) cells (Su and Gudas, 2008), 
despite its values, has been hampered by the observation of 
artifactual ligand responsiveness of the expressed RAR 
subtypes. Thus, we are presently facing a situation in which 
significant knowledge has been accumulated about the very 
early steps in retinoid action and the (patho) physiological 
impact of RAR and RXR signalling. However, what has 
remained entirely enigmatic is how a single compound upon 
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activating subtype-specific RXR-RAR heterodimers can set up 
the temporal order of complex signalling networks that are 
at the basis of (patho) physiological phenomena. 

Knowledge about the early events in retinoid signalling has 
been derived mainly from in vitro models like F9 EC cells, 
which differentiate into primary endo dermal-like cells upon 
exposure to all-trans retinoic acid (ATRA); this differentiation 
is well characterized by morphological changes and marker 
expression. F9 cells display a very low rate of spontaneous 
differentiation, such that homogeneous cell populations can 
be generated during ATRA-induced differentiation. Previous 
studies demonstrated that, while different RXR-RAR isotype 
combinations control the expression of different target genes, 
the RXRa-RARy heterodimer is essential for inducing differ- 



entiation (Taneja et al, 1996; Chiba et al 1997a, b). Together, 
these data support a model in which various RXR-RAR 
heterodimers regulate subtype-selective gene programs, of 
which RXR-RARy establishes a path that leads to the changes 
which specify a differentiated F9 cell. 

Here, we have addressed the question of how RXRa-RARy 
upon activation by ATRA sets up a sequence of temporally 
controlled events that generate different subsets of primary 
and secondarily induced gene networks. We hypothesized that 
these networks required temporally defined step(s) of diversi- 
fication, thereby forming separable gene cohorts that consti- 
tute the various facets of differentiation, such as altered 
proliferation, cell physiology, signalling, and finally terminal 
apoptogenic differentiation. To this aim, we performed RARy 



Box 1 Integrative 'omics' approach to construct the dynamic RXRa-RARy signalling network during ATRA-induced F9 
cell differentiation. 
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Cell differentiation was studied over 48 h after ATRA induction by establishing dynamic transcriptomics and ChlP-seq profilings to correlate genome-wide RXRa 
and RARy chromatin binding patterns with gene expression. RXRa and RARy metaprofiles, constructed from the cumulation of ChlP-seq patterns at all time 
points (0, 2, 6, 24, and 48 h) were instrumental for curation of the spatio-temporal binding information before integration of transcriptomics data. Combined data 
sets were used for the identification of putative RXRa-RARy target genes. In addition, the information obtained from temporal transcriptomics data sets generated 
with RAR isotype-selective agonists were incorporated in the analysis. The temporal transcription regulation information, the RXRa-RARy direct target 
annotations and presently available TF binding site annotations were integrated into the Dynamic Regulatory Events Miner (DREM) to identify decision points that 
define a co-expression regulatory map and predicted TF-based key decisions that lead to the temporal establishment of subprograms during differentiation. Finally, 
this dynamic regulatory map enabled the reconstruction of an RXRa-RARy signalling network from functional co-citations. t*h, transcriptome at time poinfh; p*h, 
chromatin binding at time-poinfh; TF, transcription factor. 
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and RXRoc chromatin immunoprecipitation (ChIP) analyses 
coupled with massive parallel sequencing (ChlP-seq) together 
with the corresponding microarray transcriptomics at five 
time points during differentiation (Box 1). To understand the 
dynamics of ATRA-regulated gene expression during differ- 
entiation, gene-regulatory decisions were inferred in silico 
from characterized targets of RXRoc— RARy and other anno- 
tated TFs (Ernst et al, 2007). This dynamic regulatory map 
was used to reconstruct RXRa-RARy signalling networks by 
integration of functional co-citation. Altogether, we present 
a genome-wide view of the temporal gene-regulatory events 
elicited by the RXRa-RARy during F9 cell differentiation. 

Results 

Genome-wide characterization of RXRa-RARy 
binding sites during ATRA-induced F9 cell 
differentiation 

We first confirmed the induction of markers {Rarfi, Hoxal, 
and Col4al) for F9 cell differentiation by RT-PCR (Supple- 
mentary Figure SI A) and the detection of binding at previ- 
ously described RAREs in the Cyp26al promoter (Loudig 
et al, 2000, 2005) using anti-RXRa antibodies (Rl and R2 in 
Supplementary Figure SIB and C). As expected, these sites 
were empty in F9 cells lacking RXRoc (Rxra _/_ ) . 

We reasoned that combining uniquely aligned reads from all 
ChlP-seq time points (0, 2, 6, 24, and 48 h) would generate a 
valuable meta binding site profile for subsequent analyses, as 
it (i) cumulates all stable and transient binding events over the 
48-h period and (ii) increases the peak calling confidence due 
to the combination of five data sets. Therefore, uniquely 
aligned reads from the RXRoc and RARy ChlP-seqs at different 
time points were combined and processed (see Materials and 
methods) to generate the corresponding metapro files . 

To identify chromatin sites occupied by RXRa-RARy hetero- 
dimers, binding sites for the two receptors in the metaprofiles 
were compared at different P- value thresholds and the percentage 
of co-occupancy was plotted for each receptor (Figure 1 A) . This 
analysis identified an optimal confidence threshold (CT40; 
P- value 10 ~ 4 ) for which all 4281 identified RARy meta sites 
were co-occupied by RXRoc. For the same CT RXRoc bound 
to 9065 additional sites, most likely as heterodimer with 
partner (s) other than RARy. Note that the implication of other 
RXRoc heterodimers in ATRA-induced F9 cell differentiation 
has been reported (Chiba et al, 1997a). 

Highly dynamic binding of RXRa-RARy during 
differentiation 

Temporal analysis of RXRoc and RARy at its 4281 meta binding 
sites revealed a highly dynamic binding (Supplementary 
Figure S2). In absence of ATRA, 2158 of the meta binding 
sites were co-occupied by RXRoc and RARy. Two hours later, 
1124 additional meta sites were occupied by the heterodimer, 
thus increasing the number of co-occupied sites; a similar 
addition of new heterodimer binding sites was observed at 
later time points, albeit with decreasing tendency (Figure IB). 
Importantly, the number of RARy-RXRoc binding sites decreased 
when cells moved through the differentiation program from 
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initially - 2000 sites at 0 h to < 1000 sites at 48 h. At 2 and 6 h, 
the gain in heterodimer binding compensated the loss of 
sites present at 0 h, while after 6 h there was an overall loss 
of RXRa-RARy binding and at 48 h only 814 were observed. 
A similar loss was observed for the number of sites that were 
newly added at a given time point and decreased thereafter. 

The observed decrease of RARy-RXRoc binding sites during 
differentiation could be due to (i) dissociation of both 
heterodimer subunits or (ii) replacement of the RXRa-RARy 
by another RXR heterodimer. Monitoring the fraction of RXRa- 
bound sites to which RARy is bound revealed that exposure to 
ATRA significantly decreased RARy co-binding to RXRa-bound 
sites over time (Figure 1C). An example is the binding of the 
RARy-RXRa heterodimer to the well-known RARE of the Rarfi 
promoter for which the level of RARy binding decreases over 
time while RXRa binding is maintained, if not increased 
(Figure ID). Most importantly, reChIP experiments, in which 
RARy or RARa is immunoprecipitated from the RXRa ChIP, 
demonstrated an unexpected strong increase of RARa co- 
occupancy at 48 h which was not observed at earlier time 
points (Figure IE and F) . Note the Rary _/ ~ and Rara _/ ~ F9 cell 
control ChlPs, which reveal the background of the assay. 

Together, the above data give not only a global view of 
the chromatin binding dynamics of the RXRa-RARy hetero- 
dimer but also provide moreover evidence for its replace- 
ment during F9 cell differentiation by RXRa heterodimers 
with other partners at common response elements. At present, 
we cannot distinguish between swapping of RXRa partners, 
i.e., dissociation followed by the formation of a distinct RXRa 
heterodimer, and the replacement of RXRa-RARy by other 
pre-formed RXRa heterodimers. 



RXRa-RARy co-occupancy correlates with gene 
induction while gene repression is largely 
independent of this heterodimer 

Transcription profiling using microarrays performed at the 
same time points as ChlP-seqs revealed a biphasic global gene 
induction with peaks at 2 and 48 h, reminiscent of results 
obtained by co-exposure to ATRA and cAMP (Harris and 
Childs, 2002). Indeed, 2h after ATRA induction 281 genes 
exhibited an induction of ^ 1.8-fold relative to Oh, followed 
by a progressive decline until 24 h (6h, 189 genes; 24 h, 128 
genes; Figure 2A). In contrast, a strong 'wave' of gene induc- 
tion was apparent at 48 h, with 926 genes getting induced. 

When comparing the differential gene expression with the 
location of RXRa or RARy inferred from the metaprofiles we 
found that >50% of the genes induced during the first 24 h 
presented an RXRa or of RXRa-RARy site within 10 kb distance 
(referred to as 'putative target genes'). Similarly as for the 
oestrogen receptor (Carroll et al, 2005, 2006), - 70% of RXRa 
(heterodimer) binding sites are beyond this distance at all time 
points and may regulate non-annotated transcripts, such as 
ncRNAs, or cognate targets through chromosomal looping 
(Supplementary Figure SID and E). At 48 h, the fraction of 
genes with RXRa/RXRa-RARy sites dropped to 34% of all 
induced genes. This reveals that the majority of gene induc- 
tions at this time are due to secondary responses. Less than 
10% of the downregulated genes presented a proximal RXRa 
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Figure 1 RXRa and RARy nuclear receptors present a highly dynamic binding to chromatin during ATRA-induced F9 differentiation. (A) Uniquely aligned reads 
sequenced from samples associated with the different time points were combined and processed to generate a meta-binding profile. The percent of RXRa and RARy 
co-occupancy relative to the total number of binding sites in their corresponding metaprofile is illustrated for different P-value confidence thresholds (CT=-10 x log 
(P-value)). The inset (Venn diagram) shows that at CT=40 all identified RARy sites are found co-occupied with RXRa. This subset of binding sites is considered 
bona fide RXRa-RARy heterodimer binding sites and has been used for all further analysis. (B) The RXRa-RARy binding sites identified in (A) are illustrated in the 
context of their temporal recruitment, duration of occupancy and dissociation (CT25). RXRa-RARy co-occupied sites per time point are subclassified based on their 
recruitment intervals and depicted by colour coding. (C) Progressive loss of RARy but not of RXRa from chromatin binding sites during ATRA-induced differentiation. 
For each time point, the fraction of RXRa-RARy co-occupied sites relative to those bound by RXRa is represented for two CT values. (D) Examples of ChlP-seq 
profiles revealing the divergent temporal binding of RXRa and RARy to the flarp promoter region; the corresponding metaprofiles (bottom panels) and the MeDiChl- 
predicted P-values (heatmaps at the right of each profile) are indicated. (E) ReChlP-qPCR quantification for temporal pattern of RXRa (primary IP) and RARy 
(secondary IP) colocalization at the Rarb promoter. flary _/ ~ cells treated with ATRA during 48 h were used to define the background. (F) ReChlP-qPCR as in (E) 
but using anti-RARa antibodies for the secondary IP; flara~ /_ cells were used as background control. In (E) and (F), the fold occupancy levels were calculated relative 
to a chromatin region localized at 18 kb downstream of Hoxbl, which corresponds to a 'cold' region. 



or RXRa-RARy binding site, suggesting that this heterodimer 
functions predominantly as positive regulator of transcription 
in this context. 
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A comparison of induced mRNA levels and gene-proximal 
temporal binding of RXRa-RARy indicated a significant corre- 
lation between binding and transcription activation. Indeed, 
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Figure 2 Temporal correlation between RXRa-RARy heterodimer binding and transcriptional regulation of putative target genes. (A) Genes exhibiting ATRA-induced 
or repressed mRNA levels at the indicated time points during F9 cell differentiation (induced genes ^ 1.8-fold; repressed genes < 0.5-fold relative to vehicle) were 
classified as putative target genes if gene-proximal RXRa or RXRa-RARy binding site was present in the CT40 metaprofiles. (B) Top panel: ranking of putative RXRa- 
RARy target genes according to the mean of their mRNA expression levels over all four time points relative to 0 h. Bottom panels: illustration of putative RXRa-RARy 
target genes ranked as above (green, relative mRNA levels) at each of the five time points during differentiation, overlaid with a display of RXRa and RARy co-binding at 
each target, expressed as the product of the corresponding confidence factors (proportional to P-value) (red for genes with fold induction levels ^ 1 .8; otherwise grey). 
(C) RNA polymerase II enrichment at TSSs and gene bodies as assessed by POLYPHEMUS from ChlP-seq assays at the indicated time points and expressed relative 
to the 0-h sample. The top 50 genes, ranked according to Poll I enrichment at their TSSs, are depicted (heatmap range ± 2a standard deviation). Note that the top 10 
genes are significantly enriched for TFs. 



sorting of putative RXRa-RARy target genes by induction 
levels revealed that at 2 h RXRa-RARy is bound predominantly 
to strongly induced genes (Figure 2B). At 6h, RXRa-RARy 
binding is more prevalent at moderately induced genes, while 
at 24 and 48 h the number of binding events in gene-proximal 
RXRa-RARy sites has dramatically decreased and the remaining 
subset is progressively associated with weakly induced genes. 

To further assess the connection between RXRa-RARy 
binding and transcription regulation of putative target genes, 
we mapped RNA Polymerase II (PolII) recruitment during 
ATRA-induced F9 cell differentiation by ChlP-seq. This 
analysis provided information about binding of PolII at both 
Transcription Start Sites (TSSs) and gene bodies (Supplemen- 
tary Figure S3). For this, the PolII binding profiles were 
processed with POLYPHEMUS (Mendoza et al, submitted), 
which entails non-linear normalization of PolII enrichment 
of multiple ChlP-seq data sets. Genes presenting proximal 
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binding sites for RXRa-RARy were subsequently ranked by 
their PolII recruitment to TSSs at a given time point relative 
to Oh. Interestingly, most of the top 50 genes (Figure 2C) 
presented significant PolII enrichment in both gene body and 
at the TSSs, indicative of active transcription. Furthermore, 
except Cyp26al and Prrl4 the top 10 genes are TFs, supporting 
a hierarchical model of ATRA-regulated gene networks in 
which RXRa-RARy induces TFs, which in turn induce their 
cognate gene programs. 



The spatio-temporal binding of RXRa and RARy 
and target gene profiling reveal distinct classes 
of temporally controlled gene induction patterns 

To link the binding of RXRoc and RARy to transcription 
activation, we clustered the putative target genes by their 
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Figure 3 Temporal (transcription) regulation defines distinct classes of RXRa-RARy target genes. (A) SOTA classification of putative RXRa-RARy target genes 
according to the indicated criteria for RXRa and RARy binding, co-binding and gene induction reveal four different classes: (i) early induced genes displaying sustained 
expression over 48 h; (ii) early but transiently induced genes; (iii) early-late transiently induced genes and (iv) late induced gene expression. Only genes that show 
coordinate heterodimer binding and gene activation at least at one time point are considered. (B) Illustration of putative target genes per class. Genes in bold were 
previously described as ATRA responsive. Heatmaps on the left (black-yellow gradient) give the P-value confidence for RXRa and RARy binding to each gene in 
the metaprofiles. Genes with more than one RXRa-RARy binding site appear several times; genes in red are validated by ChlP-qPCR and reChlP-qPCR in (D, E). 
(C) Examples of ChlP-seq profiles per class. RXRa (red) and RARy (blue) profiles are overlaid and depicted per time point. Heatmaps in the right display P-value 
confidence as in (B). (D) ChlP-qPCR validation of RXRa and RARy binding depicted as fold occupancies relative to a 'cold' region. (E) ReChlPs to assess co-binding 
of RXRa with RARy (black line) or RARa (dashed line). 



temporal receptor binding and gene expression characteristics 
using a self-organization tree algorithm (SOTA; Figure 3) . This 
classification revealed the existence of four classes of genes, 
which differ in the timing of heterodimer binding and gene 
induction (i) early induced genes with sustained expression 
over 48 h; (ii) early transiently induced genes; (iii) early-late 
transiently induced genes and (iv) late induced gene expres- 
sion (Figure 3 A and B). These classes contain several 
established RXR-RAR targets, such as Cyp26al , Rarfi or Hoxal 
(Supplementary Table I). Note that we found a third RXRa- 
RARy binding site (R3) localized ~2kb downstream of the 
Cyp26al coding region apart from the distal (R2) and proximal 
(Rl) RAREs and detected binding sites in genes shown to 
respond to ATRA but for which no RARE is described, such as 
Stra6, Stra8, Cdxl, Aqp3, Foxa2/HNF-3, and Nostrin/mDaIP2. 

For each of the four classes the timing of coordinate binding 
and gene activation was the distinctive feature, while no 
common feature could be defined for the binding of the two 
receptors before or after this phase. Indeed, the ~2.5-kb distal 
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RXRa-RARy binding site of Aquaporin [Aqp3) (Bellemere 
et al, 2008; Cao et al, 2008) was co-occupied by both receptors 
already in absence of ATRA, while binding of RARy was 
strongly reduced at 24 h and no binding of either receptor was 
apparent at 48 h (Figure 3C and D). In addition, co-activator 
components like RAC3 and p300 were recruited to this 
site at 2h and were progressively reduced at later time 
points (Supplementary Figure S4). Notably, Aqp3 expression 
increased even after receptors/co-activators disappeared from 
the locus (Figure 3B and C; Supplementary Figure S4). As for 
Apq3, RXRa-RARy occupied the putative RARE of Notch4 
(Uyttendaele et al, 1998) in absence of the cognate ligand and 
induced transcription from 2h on, but the loss of RARy 
correlated with termination of Notch4 induction and decreas- 
ing mRNA levels. In the case of Ksrl (Wang et al, 2006), 
binding of RXRa-RARy was detected at 2-6 h, followed by a 
short pulse of transcriptional induction around 6h, which 
ceased before 24 h together with the loss of receptors from 
the binding site. The late induced Nostrin (Cho et al, 1999; 
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Cho and Park, 2000) exhibited a strongly delayed binding 
of RXRa and RARy at 48 h which correlated with late RAC3 
and p300 co-activator recruitment and late gene induction. 
The RXRa-RARy co-occupancy of these binding sites at 
different time points was confirmed by reChIP assays 
(Figure 3E). In summary, the spatio-temporal cross-compar- 
ison between RXRa-RARy binding and transcriptional activa- 
tion revealed the existence of at least four different gene 
classes with distinct temporal inductions. 

The putative RXRa-RARy target genes contain a 
subset of promiscuously regulated genes that 
respond to other RAR isotypes 

To assess the selectivity and promiscuity of RAR isotype 
signalling the use of isotype-selective ligands (de Lera et al, 
2007) in the context of wild-type cells appeared to us superior 
to the use of RAR isotype-deficient cells, as such cells may 
exhibit artifactual ligand responses (Chiba et al, 1997a, b). 
To reveal RAR isoform-selective transcription of putative 
RXRa-RARy target genes, we thus used the RARy-selective 
ligand BMS961. Notably BMS961, which suffices to drive F9 
cells into differentiation (Taneja et al, 1996; see Supplementary 
Figure S5A and B) , activated 62 % of the ATRA-induced putative 
RXRa-RARy targets (Figure 4). The RARa or RAR(3-selective 
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Figure 4 RXRa-RARy putative target genes activated by specific RAR 
agonists. mRNA expression heatmaps of putative RXRa-RARy target genes 
illustrate their induction in presence of ATRA or the indicated RAR isotype- 
selective ligands. 
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BMS753 and BMS641, which do not induce F9 differentiation 
(Taneja et al, 1996 and our unpublished results), still activated 
40 and 10%, respectively, of the ATRA-induced transcrip- 
tome, thus providing evidence for both RARy selectivity and 
RAR isotype promiscuity of RXRa-RARy target genes in the 
context of F9 wild-type cells. That 38% of the ATRA-induced 
RARy-RXRa target genes were not activated by BMS961 
indicates that they are not required for F9 cell differentia- 
tion according to generally used criteria (Supplementary 
Figure S5). Mechanistically, these genes may be activated 
through direct or indirect action of RARa and/or RAR(3 
isotypes. Possible scenarios are that both RARy and RARa 
or RAR(3 heterodimers sequentially or coordinately bind to 
their regulatory regions, or that RARa or RAR(3 activate factors 
that synergize with RARy action. 



A dynamic regulatory map for ATRA-induced F9 
cell differentiation 

The above results reveal that the putative RXRa-RARy gene 
program suffices to trigger primitive endo dermal F9 cell 
differentiation. It is reasonable to assume a hierarchical 
architecture of this program in that a few key genes coordinate 
cascades of gene-regulatory events thus establishing subpro- 
gram networks. Indeed, the induction of multiple TFs supports 
a concept in which regulatory decisions are taken, albeit 
not exclusively, through TF action at defined time points. 
To identify these decisions, we used ATRA-induced temporal 
gene expression, TF-target gene annotations (NCBI database 
annotations and/or Matlnspector predictions; Cartharius et al, 
2005) and the identified putative RXRa-RARy target genes 
as input into the Dynamic Regulatory Events Miner (DREM; 
Ernst et al, 2007). DREM models bifurcation points (BPs) from 
the expression of a subset of genes that diverges from the 
co-expression pattern shared with a larger population in the 
previous time frame. In addition, DREM evaluates if a co- 
expression path is enriched for genes regulated by particular 
TFs whose action may account for, or contribute to the 
predicted bifurcation. DREM predicted six different co-expres- 
sion paths from three BPs (Figure 5A). The first BP occurs 
between 0 and 2 h and results in the establishment of three 
distinct programs generating induced (orange), constitutive 
(grey or path (iv); this class gets induced late) and repressed 
(red) cohorts. The second BP subdivides the repressed path 
between 2 and 6 h. It separates one cohort that is progressively 
induced between 24 and 48 h (path (v)) from a permanently 
repressed gene set (path (vi)). A third BP between 6 and 24 h 
derives three cohorts from the induced path; one that gets 
repressed (path (iii)) and two others that are induced with 
different kinetics and mean intensities (paths (i) and (ii)). 
To support the validity of the predicted co-expression paths, 
the three gene sets originating from the first BP were classified 
by hierarchical clustering. As shown in Figure 5B, each of 
these subsets clustered into cohorts predicted by the second 
and third BP, with the exception of related paths (i) and (ii) 
which appear as one class. 

One of the advantages of DREM is the possibility to derive 
associations between TFs and predicted BPs. In agreement 
with results described above (Figures 2A and 3), DREM 
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preferentially associates RXRoc-RARy with induced paths 
(i) and (ii). In addition, target genes of TF-like members of 
the Homeobox family (e.g., Hoxal, Hoxb2, Hoxb4, HoxbS), 
Myc, Rara, Rarb, Runxl, Jun, Foxa2, Gata4, Pbxl were also 
predicted to be enriched in these cohorts (see Supple- 



mentary Figure S6 for TF enrichment scores). Note that 
the repressed path (vi) is associated with TFs like Egrl 
(Min et al, 2008) and Sox2 (Orkin et al, 2008), which are 
involved in regulating cell proliferation and stem cell pluri- 
potency, respectively. 
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The dynamics of TF-mediated subprogramming of the 
RXRa-RARy regulon is further illustrated by the temporally 
regulated expression of TFs themselves (Figure 5C; Supple- 
mentary Figure S6B) . Indeed, with the exception of genes like 
Sox2, the majority of TFs are generally induced. Of interest 
is the biphasic response of Egrl and Myc, which together 
with Sox2, is associated with class (vi) genes. Egrl and 
Myc are induced when paths (v) and (vi) separate and get 
silenced between 6 and 24 h. This suggests that not only 
enhanced transcriptional activity but also temporally regu- 
lated expression of TFs contributes to the formation of 
temporal gene programs. 

To validate the role of DREM-predicted TFs involved in 
BPs, we performed small interference RNA (siRNA) knock- 
down assays using as readout the mRNA expression of 
differentiation markers Laminin al (Lamal), Laminin pi 
(Lambl), Laminin yl (Lamcl), type IV collagen al (Col4al); 
in addition, we monitored siRNA effects on the morpho- 
logical changes associated with differentiation (Figure 5E-G). 
We also knocked down expression of Foxal, a TF that is not 
predicted by DREM but is strongly and exclusively induced 
by ATRA and BMS961 (Figure 2C; Supplementary Figure S8; 
class I). Knockdown of Hoxb2, Hoxb5, Foxal or Foxa2 (see 
Supplementary Figure S7A for silencing efficiencies) reduced 
significantly the differentiation marker expression levels 
(Figure 5E). Notably, the expression levels of Nostrin, a late 
induced direct RXRa-RARy target, Bmp2, an established 
RA target or GAPDH were not, or only marginally affected 
(Supplementary Figure S7B). Tracking transfection with fluo- 
rescent 6-FAM revealed that transfected cells were generally 
delayed (or arrested) in differentiation, while non-transfected 
cells within the same population exhibited a differentiated 
morphology (Figure 5F). Counting of blinded samples by 
two independent persons provided a semiquantitative analy- 
sis (Figure 5G), which fully supports the notion that these 
TFs have important roles in the (temporal) regulation of 
gene networks that are at the basis of ATRA-induced cell 
differentiation. 

The dynamic map derived by DREM classified the differen- 
tially regulated genes during cell differentiation in six major 
paths, which can be distinguished by the relative enrichment 
of their components according to Gene ontology (GO) terms 
(Figure 5D; Supplementary Figure S8). Indeed, while the early 
and sustained induced paths (i) and (ii) are enriched for genes 
related to embryonic morphogenesis and actin cytoskeleton 
organization, respectively, the early temporally induced path 
(iii) is enriched for genes involved in steroid/cholesterol 



metabolic processes. The late induced path (iv) is associated 
with cell adhesion, positive regulation in response to external 
stimuli while path (v) is linked to cell-cycle regulation. 
Interestingly, the repressive path (vi) is enriched for genes 
that negatively regulate cell differentiation. 



A comprehensive ATRA-induced RXRa-RARy 
signalling network 

With the aim of enhancing the dynamic landscape of the 
RXRa-RARy regulome inferred by DREM, we reconstructed 
the corresponding gene networks on the basis of functional 
co-citation (Genomatix Bibio sphere PathwayEdition) and the 
identification of essential nodes by topology-based scoring 
methods (cytoHubba; Lin et al, 2008). The illustration of the 
resulting RXRa-RARy regulome (Figure 6; Supplementary File 
SI) depicts the relevant components of the six co-expression 
classes (compare Figure 5) and specifies their intraclass and 
interclass co-citation interactions. 

Several general features can be extracted from this dynamic 
network of co-expression classes. First, each class is unique in 
expressing a particular set of genes with similar general 
functionality, such as the TF-rich class (i). Second, genes 
regulating complex biological phenomena may appear in 
different classes with distinct expression profiles, as the 
subsequent inductions of cyclins and cyclin-dependent kinase 
inhibitors. Third, the present ChlP-seq data identify putative 
RAREs in a great number of genes, some of which are known to 
respond to retinoids (see Supplementary Table I) . Fourth, the 
described F9 RXRa-RARy regulome integrates several factors 
with important roles in other cell systems, such as Egrl and 
Notch4. Fifth, comparing regulation of the putative target 
genes by subtype-selective ligands reveals RAR subtype 
selectivity and promiscuity; moreover, the subset of genes 
commonly regulated by ATRA and BMS961 which are 
divergently regulated by RARa and RAR (3 ligands is likely 
constitute the bona fide differentiation program. 

Within class (i), topology-based scoring identified Jzin, Myc, 
Rara or Rarb as most important nodes. While the positive 
regulation of Jim and Myc expression by ATRA has been 
described (Supplementary Table I) the biphasic expression 
seen upon ATRA exposure is not maintained with the 
RARy-selective BMS961 (Supplementary Figure S6B). Indeed, 
BMS961 only recapitulates the early and late downregula- 
tion of the expression of Jun and Myc, respectively. Thus, 
the temporally regulated repression but not the induced 



Figure 5 Dynamic regulatory map of ATRA-induced transcriptome. (A) DREM co-expression analysis is represented by colour-coded paths that summarize common 
characteristics. The number of genes per co-expression path is indicated. Diamonds indicate three predicted bifurcation points (BP1-3); transcription factors (TFs) 
whose target genes are overenriched in a path are indicated. Node's size reflects the genes' expression standard deviation assigned to that node. (B) Classification of 
genes associated with the three paths generated by BP1, by hierarchical clustering of the corresponding temporal transcriptomics data leads to the subclassifications 
predicted by BP2 and BP3. (C) Transcriptional regulation of TFs associated with BP decisions. (D) Relevant Gene Ontology terms associated with each co-expression 
path. (E) mRNA expression levels of Laminin oc1 (Lamal), Laminin pi (Lambl), Laminin y1 (Lamd), type IV collagen a1 (Col4a1) in F9 cells transfected with siRNA 
constructs against TFs associated with BP3 or against Foxal , a TF induced exclusively by ATRA and BMS961 . Expression levels correspond to the mean of three 
replicates and are displayed relative to those found in GFP-control siRNA-transfected cells. (F) Morphology of siRNA-transfected cells 48 h after ATRA treatment. 
Transfected cells are identified by fluorescence from co-transfected FAM. Top panels: Hoxb2 or Foxal siRNA-transfected ATRA-treated cells. Bottom panels: mock- 
transfected vehicle-exposed undifferentiated cells and GFP siRNA-transfected ATRA-treated cells, respectively. Note that in the case of Hoxb2 or Foxal, transfected 
(fluorescent) cells are less differentiated than adjacent non-transfected cells (bar=25 |im). (G) Blinded semiquantification correlating morphological differentiation status 
and FAM-derived fluorescence by cell counting; data are the mean of two independent blinded quantifications. 
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Q ATRA responsive genes Q BMS-961/ATRA responsive genes Q TF 

^ Node's Raking score for shortest path identified using 
topology-based scoring methods (DSS) 

Figure 6 A comprehensive ATRA-RXRa/RARy signalling network. Genes associated with the different co-expression paths illustrated in Figure 5 are represented in 
the context of their functional gene co-citation interactions. For simplicity, only the top 100 hubs (coloured nodes) and their first neighbours (white nodes) are shown. 
Edge's widths correspond to the number of co-citations (limit ^5) described between nodes. Hub sizes and colours give the node's ranking based on topology scoring 
(double screening scheme of Hubba; Lin etal, 2008). This network is available in a Cytoscape format in Supplementary File S1. 



expression of these TFs correlates with cell differentiation. 
Multiple other TFs contribute to the definition of class (i). 
Apart from two other RAR isotypes, there is a strong repre- 
sentation of members of the homeobox TF family, including 
Cdxl, Meis2, some of which have well-characterized RAREs 
(Supplementary Table I) and served as validation marks 
for our ChlP-seqs. Finally, Foxal and two NR co-regulators 
[Ncoa7 and Nripl) are putative regulatory factors of class (i). 
In addition to TFs, this class contains also RA-target genes 
involved in retinoid homeostasis, including Cyp26al, Crabp2 
or Rbpl. Importantly, all of these genes are similarly regu- 
lated by ATRA and BMS961 but not by BMS753 or BMS641 
(Supplementary Figure S9), thus supporting a functional role 
in F9 cell differentiation. 

According to GO terms, class (i) is predicted to trigger 
positive regulation of transcription, cell differentiation and 
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responses to vitamin A. Class (ii), which shares a common 
ancestor with classes (i) and (hi), is characterized by the 
enrichment of genes involved in actin cytoskeleton organiza- 
tion (Supplementary Figure S8). This cohort contains also 
several apoptogenic factors, including Casp3, Casp8, Bcl2lll 
and Mcll, and the signalling factors Jak2, Rhob and Pirn; 
several of these genes are known to respond to retinoids 
(Supplementary Table I) . Comparing the induction profiles of 
these genes by the three RAR subtype-selective agonists 
indicates that their ATRA regulation may not be directly linked 
to F9 differentiation; examples for this notion are Id2, Casp 3 or 
Piml (see class (ii) in Supplementary Figure S9) . 

Several genes that are components of a similar biological 
process are found in different classes and it is tempting to 
speculate that this may be linked to their distinct temporal role 
during the differentiation process. For instance, the temporally 
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controlled expression of cell cycle-regulatory genes during 
differentiation can be rationalized from their functionalities; 
while Ccnd2 (class (ii)) expression increases, the class (hi) 
genes Ccndl, Ccnd3 and Ccnel are only transiently induced 
and regress with the expression of 'late' class (v) genes p27/ 
Cdknlb, p21/Cdknla and p57/Cdknlc, which were repressed 
at earlier time points. Note that these observations corroborate 
previous reports showing that also the protein levels of cyclin 
D2 and p27/Cdknlb increased during F9 cell differentiation, 
whereas that of cyclin Dl, D3 and cyclin E decreased (Li et al, 
2004). 

In addition to cell cycle-regulatory components, class (iii) 
contains also genes like Egrl and Notch4. The ATRA- 
responsive Egrl (Edwards et al, 1991) is sufficient to inhibit 
cell proliferation of haematopoietic stem cells (Min et al, 2008) 
where Notch4 negatively regulates cell differentiation (Vercau- 
teren and Sutherland, 2004; Ye et al, 2004). Note the temporal 
binding of RXRa-RARy to the Notch4 promoter region in F9 
cells, which correlates with its transient mRNA expression 
(Figure 3). 

The 'late' class (iv) comprises genes involved in processes 
like extracellular matrix organization [Coll8al, Sparc, Timp3) 
or cell adhesion {Cd9, Cd47, Nedd9, Cdh2, Thbsl, Sirpa, CdhS, 
Cyr61, Cdhll). These gene regulations are likely readouts of 
the morphological changes associated with differentiation, 
supported by the fact that RARa or RAR (3- specific agonists do 
not induce a similar expression as ATRA or RARy-specific 
agonists (Supplementary Figure S9). Also, late induced TFs 
like Foxa2/HNF-3, Runxl and the TF-related factor Zmizl are 
members of this class (Supplementary Table I). Finally, 
germline-specific differentiation markers, like Otx2 and Fgf4 
for ectodermal or Bmp2, Gata4 and Hnflb for endodermal 
differentiation were retrieved in classes (iv) and (vi) for which 
their temporal transcription pattern correlates with differen- 
tiation into primitive endoderm (Supplementary Figure S9). 

Taken together, the comprehensive co-expression network, 
reconstructed from an integrative analysis of RXRa-RARy 
binding and transcription regulation, illustrates the complex 
temporal coordination of a plethora of molecular processes 
during differentiation, including cell-cycle regulation, the 
activation of apoptotic/cell survival or repression of self- 
renewal/pluripotency. 



Discussion 

NR-regulated cell physiological phenomena, such as differ- 
entiation, growth, survival or death, are at the basis of 
complex physiological processes. Despite an enormous gain of 
knowledge about molecular and structural features of NR 
function (Perissi and Rosenfeld, 2005; Altucci et al, 2007; 
O'Malley and Kumar, 2009), we are still far from under- 
standing of how a single compound, such ATRA, can induce 
temporal patterns of coordinate regulation of gene networks 
that finally lead to the observed changes of cell, organ or 
animal physiology. With the availability of a number of 
genome-wide profiling technologies, it is now feasible to 
initiate integrative analyses to decipher the gene-regulatory 
phenomena and describe the dynamic regulation of gene 
networks. To this aim, we have used the ATRA-induced F9 EC 
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cell differentiation model and describe to our knowledge the 
first dynamic analysis of the RXRa-RARy regulon by an 
integrative analysis of genome-wide temporal binding of the 
RXRa-RARy heterodimer, the corresponding temporal gene 
regulation patterns and the response to subtype-selective RAR 
agonists. Our aim was to understand signal diversification and 
specification at different functional levels, ranging from the 
specific action of RAR subtypes to decisions within the process 
of differentiation that result in the formation of separate gene 
programs and finally to the description of the key components 
therein which may provide clues to the program function (s). 

One of the unexpected outcomes of this analysis is the 
astounding dynamics of RXRa-RARy binding. It was recently 
shown that ATRA-induced differentiation of mouse ES cells to 
motor neurons results in widespread changes in RAR genomic 
binding (Mahony et al, 2011) . This documented that in contrast 
to the apo oestrogen receptor, which binds to a limited number 
of target sites and acquires extensive chromatin binding 
capacity in presence of agonists (Carroll et al, 2005, 2006; 
Ceschin et al, 2011) , apoRAR already binds to a large number of 
sites (some of) which it may silence by recruiting corepressor 
complexes. On the other hand, ligand exposure was shown to 
generate de novo RAR binding to certain target genes (Chen 
et al, 1996) . Our present data not only confirm those studies at 
the genome-wide level but also demonstrate, moreover, a 
precise temporal order of gain and loss of binding sites for a 
particular RXR-RAR heterodimer during the first 48 h of 
differentiation. We observed an unexpected fluctuation of 
different RXRa heterodimers at sites initially occupied by 
RXRa-RARy. Two interesting scenarios could account for this 
phenomenon, replacement of RXRa-RARy by another RXRa 
heterodimer or swapping of RXRa partners. This temporal 
order of binding of different RAR heterodimers to the same site 
is supported by the observations (i) that reChIP experiments 
reveal a temporally distinct co-binding of RARy and RARa 
together with RXRa to the Rarb RARE (Figure IE and F) and 
(ii) that a significant number of putative RXRa-RARy target 
genes respond not only to RARy-selective but also to RAR(3 
and/or RARa-selective ligands (Figure 4; Supplementary 
Figure S9) . Taken together, our analysis of RXR-RAR binding 
reveals a highly dynamic occupancy of its target loci suggesting 
the existence of mechanism (s) that orchestrate in a clock- 
like manner the sequential recruitment and release of RXR 
heterodimers. It is tempting to speculate that site-selective 
chromatin modifications have a role in this process, which 
may involve previously observed histone methyltransferase 
and demethylase gatekeepers (Garcia-Bassets et al, 2007) . 

Such mechanism (s) may also account for the absence of 
temporal correlation between heterodimer binding and gene 
activation. Indeed, the presence of RXRa-RARy localization 
preceded in certain cases by several hours the transcription 
induction, indicating that holoRXRa-RARy binding is not 
per se sufficient for transcription induction; inversely, the 
dissociation of RXRa-RARy does not necessarily correlate with 
attenuated transcription. Due to this lack of temporal correla- 
tion the generation of metaprofiles from different time points 
facilitated the comparison with temporal gene regulation. 

The observation that distinct RAR subtype heterodimers 
can bind at identical target sites complicates the definition of 
subtype-selective target gene repertoires. As it has been shown 
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that gene deletion can generate artifactual RAR subtype 
responses (Chiba et al, 1997a, b), and affect global gene 
expression profiling even in absence of any treatment (Su and 
Gudas, 2008), we used RAR subtype-selective agonists (Gehin 
et al, 1999; Germain et al, 2004; de Lera et al, 2007) to define 
selective and promiscuous functionalities of RAR subtypes. 
Notably, the RARy-specific agonist BMS961, which is compe- 
tent to induce F9 cell differentiation, activated > 60 % of the 
putative RXRa-RARy target genes, thus confirming a causal 
link between binding and transcription activation for this 
cohort. While our analysis identified the RARy-selective gene 
program for F9 cell differentiation and revealed RAR subtype 
promiscuity, definition of the complete program will require 
the target site identification of the great number of binding 
sites that could not be linked to annotated genes. 

To approach an understanding of the decisions that lead 
to the establishment of subprograms, we used DREM for the 
in silico reconstruction of a dynamic regulatory map. This 
approach was supported by the observation that multiple TF 
were among the genes regulated by ATRA and BMS961. DREM 
analysis predicted six different gene co-expression paths, 
which result from three distinct bifurcations during the first 6 h 
of ATRA treatment. Note that this is the time span required for 
differentiation commitment (Levine et al, 1984). Some of the 
TFs associated with the DREM-predicted gene co-expression 
paths were shown here to have an important role during F9 cell 
differentiation (Figure 5E-G; Supplementary Figure S7), 
others were described previously (Gata4, Futaki et al, 2004; 
Sox2, Chew et al, 2005; Boer et al, 2007). 

To gain insight into the functionality of the co-expression 
gene cohorts and identify key players, we established a 
dynamic regulatory map and constructed a differentiation- 
related gene network from functional co -citation. While not 
exhaustive, this RXRa-RARy driven dynamic network pro- 
vides a global view on the regulatory events during ATRA- 
induced F9 cell differentiation. Indeed, genes associated with 
each class were shown to share similar general functions and 
genes of distinct classes displayed a coordinate temporal order 
of expression that is requisite for regulating complex biological 
phenomena. That an important number of RXRa-RARy target 
genes did not show up in the reconstructed network is due to 
their original identification in the present study, accounting 
for the lack of co-citation. It is worth pointing out that a 
comparison of gene regulation by subtype-selective ligands, 
as illustrated in the reconstructed network, was useful to 
validate their importance for differentiation and discover RAR 
subtype promiscuity. 

Overall, the current study illustrates for the F9 model the 
different levels of ATRA-induced signalling pathway diversifi- 
cation due to regulatory decisions at different levels and time 
points. Whereas other regulatory principles, like the chroma- 
tin modification status and co-regulator function and mod- 
ification (Ceschin et al, 2011) may be overlaid at any of these 
levels, the current gene network identifies potential nodes that 
act as key regulators of various subprograms of RA-initiated 
signal transduction during differentiation. It will be interesting 
to integrate other RXR-RAR components and compare other 
RA-regulated cell models to retrieve common and cell-specific 
features. Ultimately, it may be possible to predict critical nodes 
from applying computational models to reconstructed gene 
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networks (Fisher and Piterman, 2010) and extrapolate this 
information towards other model systems that mimic RAR- 
dysfunctional human diseases (Gronemeyer and Zelent, 2009; 
Mikesch et al, 2010), as targetable key factors for therapy. 

Materials and methods 
F9 mouse EC cell culture 

EC cells were cultured in Dulbecco's modified Eagle's medium 
supplemented with 10% fetal calf serum and 40ug/ml Gentamicin. 
Cells were seeded in gelatine-coated tissue culture plates (0.1%) and 
ATRA was added to the plates in a final concentration of 1 x 10~ 6 M 
at different time points. For assays involving RAR subtype-specific 
agonists, cells were incubated with BMS961 (RARy specific; final 
concentration 10~ 7 M), BMS753 (RARa specific; final concentration 
10" 6 M) or BMS641 (RAR|3 specific; final concentration 10" 6 M). 



ChIP and reChIP assays 

Cells were fixed with 1% para-formaldehyde (Electron Microscopy 
Sciences) for 30min at room temperature. ChIP assays were per- 
formed following standard conditions: chromatin sonication and 
immunoprecipitation in lysis buffer (50 mM Tris-Cl pH=8, 140 mM 
NaCl, 1 mM EDTA, 1 % Triton, 0.1 % Na-deoxycholate) complemented 
with protease inhibitor cocktail (Roche 11873580001); 2 x washes 
with lysis buffer; 2 x washes with lysis buffer containing 360 mM 
NaCl; 2 x washes with washing buffer (10 mM Tris-Cl pH=8, 250 mM 
LiCl, 0.5% NP-40, ImM EDTA, 0.5% Na-deoxycholate); 2 x washes 
with 1 x TE; elution at 65°C; 15 min in elution buffer (50 mM Tris-Cl 
pH=8, 10 mM EDTA, 1% SDS). RXRa and RARy have been immuno- 
precipitated with antibodies generated by immunization of rabbits 
with the following peptides: 

mRXRoc: PB105 (MDTKHFLPLDFSTQVNSSSLNSPTGRGC) , 

mRARy: PB288 (CSKPGPHPKASSEDEAPGGQGKRGQSPQPD). 

Polyclonal anti-RXRa and anti-RARy were purified from the crude 
serum by affinity chromatography. RNA Polymerase II (sc-9001 H-224), 
p300 (sc-584 N-15) and Rac3/NCoA-3 (sc-9119) antibodies were pur- 
chased from Santa Cruz biotechnology (Santa Cruz, CA). For RXRa, 
RARy, p300 and RAC3 6 x 10 6 cells were used per ChIP, whereas for 
RNA polll 2 x 10 6 cells were used. For reChlPs, at least four ChIP assays 
of 6 x 10 6 cells were used for the first IP. For reChlPs, the first antibody 
(anti-RXRa) was covalently linked to the sepharose protein A (Sigma 
P92424) using disuccinimidyl suberate (DSS). The covalently linked 
Ab beads were washed with ethanolamine (0.1 M), followed by glycin 
at pH=2.8 to remove non-covalently linked antibodies from the beads. 
Beads were washed with 50 mM sodium borate at pH=8.2 and PBS, 
and were incubated overnight at 4°C with the corresponding whole 
cell extract as in a regular ChIP assay. Following standard washing, 
elution was performed with 10 mM DTT (30 min, 37°C). The eluates 
from four ChlPs were combined, diluted at least 30 times with lysis 
buffer (containing protease inhibitors like in a regular ChIP assay), 
and incubated overnight with the second antibody (anti-RARy) and 
sepharose Protein A beads at 4°C. The subsequent steps were performed 
as for regular ChlPs. 

The immunoprecipitated chromatin was validated using positive 
(recruitment to known targets) and negative ('cold' region) controls 
and the binding was expressed as enrichment relative to the whole 
cell extract input control (% input) and/or relative to a 'cold' region 
(fold occupancy); validation assays were performed by quantita- 
tive real-time PCR (qPCR, Roche LC480 light cycler device) using 
Quantitect (Qiagen). 



Massive parallel sequencing 

After qPCR validation, immunoprecipitated chromatin was quantified 
using Qubit (Quant-It dsDNA HS Assay Kit; Invitrogen) . In all, 10 ng of 
the ChlPed material was used for preparing the sequencing library 
(ChlP-seq DNA sample preparation kit, Illumina) . In all, 5 pmol of the 
library was used per flow cell in the Solexa 2G Genome analyzer 
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(Illumina). The Illumina Pipeline vl.4.0 was used for primary data 
analysis (image processing: Firecrest; base calling: Bustard; alignment: 
Gerald) from which uniquely aligned reads with up to two mismatches 
relative to the mouse mm9 reference genome were kept. 



Peak detection approach 

MACS v. 1.3.6 (Zhang et al, 2008) was used as peak caller at the first 
level of data treatment to obtain signal intensity wiggle files and 
annotated peak regions using a Poisson model distribution (P-value 
confidence cutoff: 1 x 10~ 5 ). Considering that a certain number of false 
positive regions were observed during this treatment, which in our 
hands cannot be decreases without compromising true positives by 
playing with the P-value confidence parameter, we used subsequently 
an approach based on a machine-learning algorithm able to define the 
model distribution associated with the ChlP-seq data set under study. 
This approach, implemented in the R package MeDiChI, has been 
previously described (Reiss et al, 2008) for analysis of ChlP-chip 
assays; its implementation for ChlP-seq assays will be described 
elsewhere (manuscript in preparation) . 



RT-PCR and microarrays 

Total RNA was extracted from cells, treated with ATRA during different 
period of times, using the GENEElute™ RNA extraction kit (Sigma) . 
In all, 2(ig of the eluted RNA was used for reverse transcription 
(AMV-RTase, Roche; oligo(dT) New England Biolabs; 1 h; 42°C). 
The cDNA was diluted 10-fold and used for real-time qPCR (Roche 
LC480). Expression of the following marker genes was assessed to 
follow the process of F9 cell differentiation: 

For early response to the differentiation inductor (ATRA) : 

Hoxal (CCC AGACGGCTACTTACC AG ; CATGGGAGTCGAGAGGTTTC) ; 

Rarb (GATCCTGGATTTCTACACCG; CACTGACGCCATAGTGGTA) . 

For late response to the differentiation inductor (ATRA) : 

CoUal (ATGCCCTTTCTCTTCTGCAA; ATCCACAGTGAGGACCAACC); 

Lamal (CCGACAACCTCCTCTTCTACC; TCTCCACTGCGAGAAAGTCA); 

Lambl (TCTATGCTCGGCAGTGTGAC; CAGTGGTCTCCTGACCCAAT); 

Lamcl (GGCCGAGTGCCTACAACTT; CAGTGGCAGTTACCCATTCC) . 

During data analysis, all qPCR values were normalized relative 
to the constitutively expressed 36B4 mRNA (AATCTCCAGAGGCACC 
ATTG; CCGATCTGCAGACACACACT). 

Using 250 ng of starting total RNA, biotin-labelled cDNAs were 
synthesized and hybridized on Affymetrix GeneChip® Mouse Gene 
1.0 ST Array (Affymetrix, Santa Clara, CA, USA) according to the 
manufacturer's recommendations as described in 'GeneChip(g) whole 
transcript sense target labelling assay manual' (P/N 701880 Rev.4). 
The arrays were further washed and stained with streptavidine- 
phycoerythrin in an Affymetrix GeneChip(g) Fluidics station 450 
using the script protocol F450-0007, then scanned with an Affymetrix 
Gene Chip Scanner 3000 7G. Expression values were generated with 
Affymetrix software Expression Console version 1.1, using sketch 
quantile normalization and median polish summarization as in 
Robust Multiarray Analysis. 



RNA polymerase II ChlP-seq data processing 

RNA Polymerase II genome-wide intensity profiles were generated 
using MACS v. 1.3.6 (Zhang et al, 2008) as peak caller and signal 
intensity profiles were processed by MeDiChI to identify chromatin 
regions presenting significant RNA Polymerase II enrichment (P-value 
cutoff: lxlO -2 ). MeDiChl-predicted enriched regions were further 
filtered according to their genomic localization by comparing them 
with the mm9 coding region Ref-seq annotation. Thereby coding 
regions presenting significant levels of RNA PolII at the TSS were 
identified and selected for further analysis. A Quantile normalization 
procedure was used to enable a comparison between different 
samples. Normalized profiles were compared in the context of their 
signal intensities to identify local changes in the Pol II occupancy 
relative to the control sample. Note that this approach differs from 
previously described linear corrections between samples based on 
total numbers of reads; it has been implemented in the R package 



POLYPHEMUS (Mendoza et al, submitted). The final comparisons 
between the PolII levels at the TSS relative to the average behaviour 
through the gene body have been used as readouts to identify and 
describe genes presenting an enhancement of their transcriptional 
activity at the different analysed time points (Figure 2). 



Data integration 

The global RXRa and RARy localization at different time points during 
ATRA-induced differentiation has been followed by ChlP-seq. To 
increase the confidence in binding localization assessment, we 
collected the mappable reads from all five time points in a single file 
called metaprofile. The metaprofile associated with the localization of 
RXRa and RARy was processed by MACS and MeDiChI to identify 
significant enriched regions. Interestingly, this approach generates a 
higher signal/background contrast, thus increasing the sensitivity of 
peak calling. Taking in consideration the heterodimeric nature 
between RXRa and RARy, we compared the two metaprofiles to 
identify the optimal confidence parameter at which the coexistence 
between RXRa and RARy is maximal (see Figure 1A). Indeed, for a P- 
value threshold of 1 x 10" 4 (CT40; where CT=-10 x log (P-value)) all 
identified RARy sites in the metaprofile colocalize with RXRa. These 
co-occupied sites were further evaluated in the context of their 
temporal RXRa and RARy co-occupancy. Whereas a CT40 has been 
used for metaprofiles binding sites selection, the same parameters are 
too stringent for the analysis of time-point profiles. We have taken a 
CT25 per time-point profile for evaluating the RXRa and/or RARy 
occupancy at the pre-identified sites (Figure IB). 

RXRa-RARy co-occupied sites [metaprofiles comparison) were 
annotated based on their proximity (<10kb) to transcriptionally 
active coding regions (upregulation and downregulation levels relative 
to control sample: ratio cutoff ^1.8 and ^0.5, respectively); 
annotated genes are referred to as 'putative' target genes. In order to 
gain temporal information about the correlation between RXRa-RARy 
binding and putative target gene transcription, the RXRa and RARy 
binding at a given time point was compared with the temporal mRNA 
levels of putative target genes. To further remove possible false posi- 
tive annotations, genes presenting coexistence between RXRa-RARy 
binding and differential mRNA expression levels at least at one time 
point were retained during selection. For classification purposes, the 
coexistence of several events was scored per time point following a 
hierarchical order of importance in the context of gene regulation: 



Evaluated event Score 



— 0 

RXRa 1 

RARy 2 

RXRa; RARy 3 

Gene induction 4 

RXRa; gene induction 5 

RARy; gene induction 6 

RXRa; RARy; gene induction 7 



Finally, the temporal incidence of the evaluated events was 
classified using an SOTA approach (Euclidean distance, Max. 
cycles=7, cell variability P-value=0.01) under the open access 
multiExperiment Viewer (Saeed et al, 2003, 2006; see Figure 3). 



Dynamic regulatory maps and Network 
reconstruction 

RXRa-RARy putative target gene information as well as further TF- 
gene regulatory interaction annotations extracted from the NCBI 
database and/or predicted by Matlnspector (Cartharius et al, 2005) 
were combined into a binary matrix ('1' for TF-target gene association; 
otherwise '0'). Furthermore, time-point transcriptomics data were 
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expressed in log2 ratios relative to the 0-h control. These two data 
sets were integrated into the DREM (Ernst et al, 2007) to establish 
dynamic regulatory maps. In addition to classifying temporal sets of 
gene expression data into co-expression paths, DREM predicts BPs, 
defined as a time point at which the co-expression behaviour of a 
subset of genes diverges from the common path present in a previ- 
ous time point. Moreover, the association of a given TF with the 
corresponding BPs is estimated by evaluating the enrichment of its 
target genes using hypergeometric distributions relative to the genes 
associated with the common path. 

To construct gene networks, genes associated with each predicted 
co-expression path were evaluated in the context of their functional co- 
citation relationships (Genomatix Bibiosphere PathwayEdition) . To 
identify functional co-citation interactions between different co- 
expression paths, all genes used for the initial analysis were evaluated 
for their functional co-citations relationships as explained before. 
Finally, the intra- and inter-co-expression paths co-citation interac- 
tions were integrated into Cytoscape (Shannon et al, 2003; Cline et al, 
2007) . To further increase the confidence of the reconstructed network, 
topological information concerning the number of edges per node 
(Hubba; Lin et al, 2008) were used as filter in the final representation. 
Briefly, we have used Hubba's Double screening scheme (DSS) 
approach, which combines a Maximum Neighbourhood Component 
(MNC) with a Density of Maximum Neighbourhood Component 
(DMNC) to score nodes based on the number of their interconnections. 
The final gene-network representation corresponds to the top 100 
ranked nodes (red to green colour code) as well as their corresponding 
first level neighbours. Additional information, like TF annotations, 
BMS961 and ATRA responsiveness, number of functional co-citations 
between nodes (edge's broadness) and co-expression path belonging 
(colour coded) are also included in the final display of this RXRa-RARy 
induced genes network, which is also available in Cytoscape format 
(Supplementary File SI). Transcriptomics data associated with ATRA 
and RAR subtype-specific agonist treatments were also included as 
attributes for the reconstructed network (Supplementary Files S2-S6). 



siRNA transfections 

F9 cells were transfected with siRNA oligomers directed against 
RNA target sequences for the following TFs: Hoxb2 (QIAGEN; 
FlexiTube GeneSolution: Cat. No. SI01069089, SI01069096, SI01069075, 
SI01069082; working concentration: 50 nM); Hoxb5 (QIAGEN; FlexiTube 
GeneSolution: Cat. No. SI01069173, SI01069180, SI01069159, SI01069166; 
working concentration: 50 nM); Foxal (QIAGEN; FlexiTube GeneSolution: 
Cat. No. SI01004493, SI01004500, SI01004479, SI01004486; working 
concentration: 50 nM); Foxa2 (QIAGEN; FlexiTube GeneSolution: 
Cat. No. SI01004528, SI02737182, SI01004514, SI01004521; working 
concentration: 50 nM); Gata4 (QIAGEN; FlexiTube GeneSolution: 
Cat. No. SI01009813, SI01009820, SI01009799, SI01009806; working 
concentration: 50 nM). In addition, F9 cells were transfected with an 
siRNA oligomer directed against the RNA target sequence of GFP 
(Dharmacon; P-002048-01-20; working concentration: 50 nM) as a 
control for the specificity of the assay. The transfection efficiency has 
been followed by co-transfecting the previous oligomers with the 
6-FAM labelled siGLO transfection indicator (Dharmacon; D-001630- 
01-05; working concentration: 50 nM). F9 cells were transfected with 
Lipofectamine 2000 (Invitrogen; Cat. No. 11668-027) during 18 h, then 
medium has been changed and cells were treated either with ATRA 
or with ethanol during 48 h as previously described. 



Data access 

Affymetrix microarrays as well as Illumina platform ChlP-seq data 
described in this study are available from the Gene Expression 
Omnibus database (http://www.ncbi.nlm.nih.gov/geo) under accession 
number GSE30539. 



Supplementary information 

Supplementary information is available at the Molecular Systems 
Biology website (www.nature.com/msb). 
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